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Abstract 

We focus on the Gibbs free energy AG for nucleating a droplet of the stable phase (e.g. solid) 
inside the metastable parent phase (e.g. liquid), close to the first-order transition temperature. This 
quantity is central to the theory of homogeneous nucleation, since it superintends the nucleation 
rate. We recently introduced a field theory describing the dependence of AG on the droplet 
volume V, taking into account besides the microscopic fuzziness of the droplet-parent interface, 
also small fluctuations around the spherical shape whose effect, assuming isotropy, was found to be 
a characteristic logarithmic term. Here we extend this theory, introducing the effect of anisotropy 
in the surface tension, and show that in the limit of strong anisotropy AG (V) once more develops 
a term logarithmic on V, now with a prefactor of opposite sign with respect to the isotropic case. 
Based on this result, we argue that the geometrical shape that large solid nuclei mostly prefer could 
be inferred from the prefactor of the logarithmic term in the droplet free energy, as determined 
from the optimization of its near-coexistence profile. 
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I. INTRODUCTION 



When a homogeneous, defect-free bulk system is brought across a first-order phase bound- 
ary, it may survive in its metastable state even for a long time, until the stable phase sponta- 
neously nucleates [HE]- The nucleation process has attracted much attention over the years, 
both from a fundamental point of view as well as for its great practical interest. To mention 
but one example, a better control of crystal nucleation in protein solutions could help hin- 
der protein condensation which is at the heart of several human pathologies [3]. Thermal 
fluctuations continuously sprout droplets of the stable phase inside the metastable mother 
phase. Small droplets dissolve, for the gain in volume free energy fails to compensate the 
loss in surface free energy. Occasionally, a droplet is sufficiently large that it is favorable for 
it to grow. Once this happens, the solid nucleus expands until the whole liquid crystallizes. 
Quenching the system deeper and deeper lowers the nucleation barrier until the point where 
the barrier vanishes (kinetic spinodal limit). Beyond this threshold, nucleation ceases and 
the phase transition occurs through spinodal decomposition and coarsening (i.e., uniformly 
throughout the material). Classical nucleation theory (CNT) jlHO] provides the simplest the- 
oretical framework in which the initial stage of the phase transformation can be described. 
In this theory, an isolated droplet is schematized, regardless of its size, as a sphere of bulk 
solid, separated from the liquid by a sharp interface with a constant free-energy cost per 
unit area a ("capillarity approximation"). This gives rise to a (Gibbs) free-energy difference 
between the supercooled liquid system with and without a solid cluster, that is 

AG(V) = -p s \Afi\V + (367i) 1/3 aV 2/:i , (1.1) 

where V is the cluster volume, A/x < is the difference in chemical potential between solid 
and liquid, and p s is the bulk-solid number density. The droplet grows if it exceeds a critical 
size V* corresponding to the maximum AG (= AG*), which thus provides the activation 
barrier to nucleation [7]. 

The cluster free energy AG(V) can be accessed numerically via the statistics of cluster 



size, through which the validity of Eq. (1.1) for specific model interactions can be directly 
tested. We recently showed [8j that the accuracy of CNT is less than satisfactory in estimat- 
ing the size probability distribution of clusters, especially the smaller ones, implying that 
interface-tension estimates based on the use of CNT are systematically in error. We then 
proposed a more detailed field theory of the nucleation barrier, based on the assumption 



that clusters are soft and not sharp, and can deviate mildly from the spherical shape ("qua- 
sispherical" approximation). If the solid-liquid interface tension is taken to be isotropic, the 
volume dependence of the Gibbs free energy of a cluster is of the Dillmann-Meier form [9], 

AG(V) = -pJAdV + AV 2/3 + BV 1/3 + C- -k B T\n ^ , (1.2) 

9 a d 

where A, B, and C can all be expressed as explicit functions of the "microscopic" parameters 
entering a Landau free energy, and a is a microscopic length. It turned out that the numerical 
profiles of AG in a few test cases and at various supersaturations are better reproduced by 
this theory. 

Here we critically reconsider the most severe assumption made in that derivation, namely 
the isotropy of the solid-liquid interface tension. We show that the theory introduced in 
Ref. [8] can be extended relaxing this important approximation, and that the results change. 
Starting once again from a Landau-like theory, we derive an interface Hamiltonian, that 
allows estimating the probability of observing a cluster of any shape and size. The angular 
dependence of the interface tension is taken into account by terms that depend on the local 
orientation of the cluster surface. Within this framework, we calculate AG(V) in the limit 
of strong surface anisotropy and compare it with the isotropic case. For large anisotropy, 
the cluster free energy still retains at large size a logarithmic term, however with a prefactor 
of opposite sign to the isotropic one. On account of this, we suggest that the nominal 
shape of large solid nuclei could be guessed from the optimization of the actual AG(V) 
close to coexistence. Looking for a numerical exemplification, we conducted 3D Monte 
Carlo simulations of the Ising model extracting AG(V) for clusters of variable size V, at 
various distances from coexistence. Although we could not really attain sizes where the 
anisotropic shape effects are heavy, we do detect evidence that the sign is as expected for 
large anisotropy. 

The paper is organized as follows. We start in Section II by relaxing the approximation 
of an infinitely sharp cluster interface, with the introduction of a Landau free energy. From 
that, an effective sharp-interface Hamiltonian is derived in Section III, as an intermediate 
step to building up a field theory for isotropic surfaces where small shape fluctuations are 
allowed (Section IV. A). Eventually, this leads to a modified-CNT expression of AG(V). In 
Section IV. B, the issue of interface anisotropics is addressed, and we show by examples how 
the dependence of the interface free energy on the local surface normal affects the formation 
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energy of a large cluster. Next, in Section V, we check our theory against old and fresh 
Monte Carlo simulation data for the nucleation barrier to magnetization reversal in the 3D 
Ising model above the roughening temperature. While confirming that CNT is not generally 
adequate to fit the numerical AG(V) data, this analysis also gives a quantitative measure of 
the errors made with CNT and demonstrates their cancellation in the more general theory. 
Finally, our conclusions are presented in Section VI. 



II. DIFFUSE INTERFACE: LANDAU THEORY 



The main assumption behind CNT is that of a sharp and spherical cluster surface. A 
way to relax this approximation is through the introduction of a scalar, non-conserved 
order-parameter (OP) field 0(x) ( "crystallinity" ) which varies smoothly from one phase to 
the other. Hence, the solid-liquid interface becomes diffuse in space, even though only on a 
microscopic scale. In practice, <fi may be thought of as the local value of the main Fourier 
coefficients of the crystal-periodic one-body density n(x), i.e., those relative to the reciprocal- 
lattice vectors which are closest in modulus to the point where the liquid structure factor 
reaches its maximum pU] . Otherwise, may be identified with the parameter discriminating 
between solid and liquid in an ansatz like 

\ 3/2 

) 5>-^ x - R ) 2 = p s J2e- G2m) e* G -\ (2.1) 

• % ' R G 

assuming a specific crystal symmetry and an overall number density p s . 

Across the solid-liquid interface, <fi is no longer constant and, for a system with short-range 
forces, the thermodynamic cost of the interface may be described through the free-energy 
functional [TTHT51 



nix 



0fofi] = / d 3 x {^(V0) 2 + ^(V 2 0) 2 + ^(0(x))} , (2.2) 

where c, k > are stiffness parameters dependent on the interface orientation as defined by 
the unit normal n and g(<p) is the specific Landau free energy of the homogeneous system, 



taken the bulk liquid as a reference. In Eq. (2.2), besides the customary square-gradient 



term, also a square-laplacian term appears. This is the next-to-leading isotropic term in the 
gradient expansion of the Landau free-energy density [16J. Even though being a fourth-order 
gradient term, it is however only second-order in the order parameter, and this places it on 



the same footing as the square-gradient term (hence, potentially relevant). We shall see 
below that, without such a term, the bending rigidity (i.e., the coefficient of H 2 in Eq. (3.15) 



below) would simply be zero. Below the melting temperature T m , g shows, besides the liquid 
minimum, also a second and deeper solid minimum. Exactly at coexistence, the two minima 
are equal, falling at = <p S Q in the bulk solid and at <p + = in the bulk liquid, which 
means that g{4> s o) — <?(0) = while g(<ft) > otherwise. 

When boundary conditions are applied such that <p — > <fi± for z — > ±oo, a planar interface 
orthogonal to z is forced to appear in the system. The corresponding OP profile is the 



stationary solution 4>o(z; n) of (2.2) that satisfies the boundary conditions: 



c(n)$ - «(fi)C = ||(0 O ; T = T m ), with o (-oo) = sO and o (+oo) = . (2.3) 
From now on, we simplify the notation by dropping any reference to n in c, k, and </>q. 



Equation (2.3) can be simplified by multiplying both sides by 4>'q(z) and integrating by 



parts. We thus arrive at a new boundary value problem: 



K 0o0o" = ^00 + ^4>'q - 9{<Po) , with (f) (-oo) = <f) s0 and o (+°o) = . (2.4) 

Obviously, Q[4>o\ represents the free-energy cost of the interface at T = T m . 

At temperature below coexistence, the absolute minimum of g(<ft) falls at <fi = <ft s > for 
AT = T — T m < 0. This can be described by 

9{<f>) = c 2 2 + c 3 3 + c 4 4 + . . . (2.5) 

with c 2 = c 2 o + c' 20 AT (c 2 o, c' 20 > 0), all other c n coefficients being constant. 

For the remaining part of this Section, we will assume that c and k do not depend on n. 

Under this condition, a large solid cluster can be assumed to be spherical, with a OP profile 

described by <po(r — R) [13], provided the center of (fio(z) is at z = 0. From this ansatz, in 

[8] we derived an expression for the cluster free energy, 

/ 9r) L f L \ 4 
AG(R) = A7cR 2 a L — + -7rR 3 Ps \A(,\ , (2.6) 

in terms of quantities (cr L , 5 L , e L ) which depend linearly on the supersaturation |A/i| oc |AT|. 



Equation (2.6) resembles the CNT expression, Eq. (1.1), with the crucial difference that the 



interface free energy is now a function of both R and T: 

~R~ + R 2 



a(R;T) = a L (l-^- + -) . (2.7) 



Exactly of this form is the tension of the equilibrium interface between a liquid droplet and 
the vapour background in the Lennard- Jones model, as being extracted from the particle- 
number histogram in grand-canonical simulations of samples of increasing size [TTJ . At 
coexistence, the solid-liquid interface tension and the Tolman length [18] are given by: 

/+oo 
dz[c<f)' 2 {z) + 2^ 2 {z)} ; 

f +0O dzz[c(P' 2 (z) + 2K(f)'' 2 (z)] 

5 m = 5 (T m ) = -±=%- (2.8) 

j+-dzW 2 {z) + 2^' 2 {z)} 

A nonzero S m occurs if and when 4>q(z) is asymmetric around zero, as is generally the case for 



the interface between phases of a different nature (see Appendix A). Summing up, Eq. (2.6) 
describes the corrections to CNT which arise by replacing the assumption of a sharp solid- 
liquid interface with a more realistic finite width, in the case of isotropic surface tension and 
Tolman length. 

III. SHAPE FLUCTUATIONS: THE INTERFACE HAMILTONIAN 

A real cluster may be spherical only on average. Far from being static, clusters fluctuate 
widely away from their mean shape [T^l 20J. To describe fluctuations, we switch from a 
description in terms of the crystallinity OP to another in which the cluster shape itself 
rises to the role of fundamental variable. We begin by deriving a coarse-grained, purely 
geometrical Hamiltonian for the cluster surface directly from the microscopic free-energy 



functional (2.2), under the assumption of small deviations of the interface from planarity. 
The outcome is a Canham-Helfrich (CH) Hamiltonian [2D [22], containing spontaneous- 
curvature and bending penalty terms in addition to interface tension. 

For the present derivation, we build on Refs. [231 EI]- Other attempts to derive an ef- 
fective interface Hamiltonian from a mean-field density functional are described in [23 EE] ■ 
Let the cluster "surface" be depicted as a closed mathematical surface £ embedded in 
three-dimensional space and let R(u, v ) be the parametrization (coordinate patch) of an 
infinitesimal piece of E. We switch from 3D cartesian coordinates, r = (x,y,z), to new 
coordinates q a = (u, v, () (tangential and normal to E) by the transformation 

r = R(m, v ) + (ri(u, v ) , (3.1) 
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where 



R„ A R, 



n(u,v) = |p p. (3.2) 
\K U A R„| 

is the unit normal to S. For a patch that deviates only slightly from planarity, we may 
adopt a free energy Q[(f>o(((x,y, z))], thus arriving at the surface Hamiltonian 



with J 



■ (r„ A r 



du dv dc J [l (V0 o (C)) 2 + \ (v 2 o (C)) 2 + g(MO) 



(3.3) 



|n • (r u Ar„)|. In order to make Eq. (3.3) simpler, it is convenient 



to view the patch as parametrized in terms of orthonormal, arc- length coordinates, i.e., 
R u ■ R v = and |R n | = |Ru| = 1 all over the patch. Although this construction is rigorously 
possible only for surfaces having zero Gaussian curvature (K = 0) |27J, we can reasonably 
expect that only small errors of order K are made for quasiplanar interfaces. With this 
caution in mind, we go on to get (see Appendix B): 



dr 

du 
dr 

dv 
Or 



-(r g R u + (1 - C4 2) )R, 



n 



(3.4) 



where k£ ^ and Kn' are the normal curvatures of the u- and w-lines respectively, and r g = 
Tg = —Tg 2 ^ is the geodetic torsion. From Eqs. (3.4), we readily derive the metric tensor g a p, 



.(2) 



9af3 



dr dr 



dq a dq p 
and the Jacobian, 



-2(T g + C 2 T g (k 



\ 2 



(1) r «£> 



-2(T g + (\ (k$P + K 

(l - C« 



(2) 





2^2 
9 







(3.5) 



J = (i - c«L 1} ) (i - c4 2) ) - c 



\ 2 



(3.6) 



g being the determinant of (3.5). Considering that covariant and contravariant components 
of a vector are built by projecting it on the bases Vg Q and dr/dq a , respectively, we can 
calculate the gradient of a scalar field and the divergence of a vector field A in local 
coordinates as follows: 



V0 



9 



a/3 



dr 



dq a dq p 



and V • A 



J__d_ 
y/gdqc 



(3.7) 
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gap being the inverse of (3.5). In particular 



V0(C) = 0'(C)n and V 2 0(C) = <P"(Q + 0'(C)V • n 



(3.8) 



where 



V ■ n 



1 



Finally, the mean and Gaussian curvatures of the patch are given by 



H 



1, 



-V • n 



C=o 



2 { K n^ + K n ^) 



(3.9) 



(3.10) 



and 



n 



<9n <9n\ 

9m <9t> / 



K (1) K (2) 



(3.11) 

Hence, 1) the mean curvature, which is defined only up to a sign depending on our conven- 
tion on the orientation of n, is half the sum of the two normal curvatures relative to any 
orthogonal parametrization, i.e., not necessarily the two principal curvatures; 2) since K 
is the product of the two principal curvatures, the geodetic torsion must vanish when the 
coordinate lines are also lines of curvature. 



— K 



We are now in a position to simplify Eq. (3.3). Upon using Eq. (2.4) to eliminate g(4>o) m 
favor of (c/2)0' o 2 + (3k/2, v >„ 
we eventually get 



and inserting Eqs. (3.6), (3.8), (3.10), and (3.11), 



H s = j du dv dC (1 + 2(H + ( 2 K) |c0' o 2 (C) + \n<Pl 2 



(0 



+ 



fflC) + 0'o(C) ^h + Ik ) -«(^(O0g(C))' 



(3.12) 



We now argue that, to a first approximation, any term of order higher than H 2 and K can 
be discarded. Moreover, J dudv = J dS since |R U A~R V \ = 1. Lastly, the geodetic torsion 
vanishes if we perform a change of integration variables (that is, a change of parametrization) 
such that the coordinate lines are also lines of curvature [28]. In the end, we are left with 
the classic Canham-Helfrich Hamiltonian for fluid membranes: 



U s = J dS (a + bH + cH 2 + dK) 



(3.13) 
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with the following explicit expressions for the coefficients: 

dcK 2 (O + 2/#' /2 (C)] ; 
dCCK 2 (C) + 2/#'o /2 (C)] ; 
dC0o 2 (C); 

-oo 

o = I dc {C 2 K 2 (0 + 2«^g a (0] - < 2 (0} 



+00 



—00 
+00 



2k 



(3.14) 

A few remarks are now in order: 1) H and K are reparametrization invariants, hence no 
ambiguity arises from the arbitrariness of the parametrization used. 2) The above derivation 
actually applies for just one E patch. However, upon viewing X as the union of many 



disjoint patches, the Hamiltonian (3.13) holds for the whole £ as well. 3) As anticipated, 



the coefficient d of the K term in (3.13) could be different from the quoted one since a 



parametrization in terms of orthonormal coordinates does not generally exist. However, as 
far as we only allow for clusters with the topology of a sphere, J s dS K takes the constant 
value of Air by the Gauss-Bonnet theorem and the K term in l-i s can be dropped. Upon 



comparing the definition of a and b in Eqs. (3.14) with Eqs. (2.8), we can rewrite Eq. (3.13) 



in the form (restoring everywhere the dependence upon interface orientation): 



dS (a m (h) - 2a m (h)5 m (h)H + 2X(h)H 2 ) 



(3.15) 



where A = c/2 (we note that A = K<fi 2 s0 / (?>€) under the same hypotheses for which Eq. (A. 15) 
holds). 4) The term linear in H is related to the spontaneous curvature of £, Hq = — 6/ (2c), 
which is proportional to the Tolman length 8 m . A nonzero value of H yields a difference 
in energy between inward and outward interface protrusions, thus entailing a non-zero 5 m . 
The additional fact that in systems, such as the Ising model, where the symmetry is perfect 
between the two phases then 5 m = 0, has long been known [T3] . 



We point out that Eq. (3.15) retains the same form as in the isotropic case [8j. In the 



general anisotropic case, the dependence of the Hamiltonian parameters on the interface 
normal is through the constants c and K, and the function <f>o(z). 
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IV. THE CLUSTER FREE ENERGY IN TWO EXTREME CASES: 
ISOTROPIC AND STRONGLY ANISOTROPIC INTERFACE TENSION 



Considering that every single realization of the profile of the cluster surface should be 
sampled in equilibrium with a weight proportional to exp{— 0% s }, it is natural to define a 
volume-dependent cost of cluster formation through 

AG(V) = - Ps \Afi\V + F S (V) (4.1) 

with 



F a (V) = -k B T\nZ s (V) = -k B T\n^a 3 J £>£ e~ m ° 6{V[E] - V) \ . 



(4.2) 



In the above expression of the constrained partition function Z s , a — p s 1//3 is a microscopic 
length of the system, V[S] is the volume enclosed by the closed surface E, and VT, a yet-to- 
be-specified integral measure. 



While the calculation of F s for a realistic form of n-dependent parameters in (3.15) is 
certainly possible numerically once the admissible surfaces have been parametrized in terms 
of a basis of eigenfunctions, some restrictions are to be made in practice if we want to make 
analytical progress. In the following, we examine two limiting cases for <x m (n), according 
to whether it is constant or strongly anisotropic. In general, a strongly anisotropic <x m is 
typical of e.g. systems where melting is very strongly first order, implying very sharp and 
thus direction dependent solid-liquid interfaces, such as for example in the case of alkali 
halides [29J. That brings about a non-spherical cluster shape through the prescription that 
the surface free energy be the minimum possible for the given cluster volume V. The same 
condition is responsible for a spherical shape when the interface free energy is isotropic. 



A. Isotropic interfaces 



If (J m ,S m , and A in Eq. (3.15) do not depend on n, the shape of a cluster is on average 



spherical. We here compute the free energy (4.2) assuming small deviations from this shape. 

Neglecting overhangs and liquid inclusions, let r = R(6, 0) be the equation of X in 
spherical coordinates. We assume only small deviations from a sphere, i.e., R(6, <f)) = Rq[1 + 
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e(0, </>)], with e(8,4>) -C 1 [30]. Then, we expand e(6,4>) in real spherical harmonics, 

oo I 

<°^) = J2 £*i,m*WM), (4-3) 

1=1 m=—l 

and we agree to ignore, from now on, all terms beyond second-order in the coefficients 
Xi,m- With these specifications, we obtain approximate expressions for the area of £ and its 
enclosed volume, as well as for the mean curvature H. Upon inserting this form of 1-L S in 



terms of the x^ m into Eq. (4.2), we are left with the evaluation of a Gaussian integral. While 
we refer the reader to Appendix C for all the technicalities, we here quote the result of the 
calculation. The free energy cost of cluster formation for large V is 



AG(V) = -p s |A yU |\/ + (367r) 1 /V Q V 2 / 3 -(3847r 2 ) 1 / 3 ( T QS ( 5 Q V 1 



/3 



r QS e QS_ lk B T\n ^vr) 1 / 3 (- 



+ Ana^e^ - -k B T\n | (36n)^ ( ' ) , (4.4) 



where cr^ 3 ,^ 3 , and e^ s can be read in Eq. (C.21). The above formula is strictly valid 



only near coexistence, where the various assumptions beneath its derivation are expected to 
hold true. We have thus found that the surface free energy has a form consistent with the 
Dillmann-Meier ansatz, with T-dependent parameters a^jS^, and that are different 
(even at T m !) from the corresponding ones in Landau-theory a L ,5 L , and e L , and with a 
universal logarithmic correction to the mean-field form of AG. This term is responsible for 
the well known i?* 7//3 exponential prefactor to the nucleation rate |31j . 



B. Anisotropic interfaces 



We now consider an interface tension of the form: 



o n 



°"ioo 



(4.5) 



with M oo, written in terms of the cartesian components of the outer normal to the 
cluster surface. In the infinite-M limit, the equilibrium crystal shape is a cube, though 
rectangular cuboids are also admissible, though not optimal, shapes (they arise at non-zero 



temperatures). The terms in Eq. (3.15) beyond the first are singular in the M — > oo limit; 



however, they would contribute to the surface free energy if M were large but not infinite, 
see more in Appendix D. In the same Appendix we show that the asymptotic, large-V 
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free-energy cost of cluster formation is given by: 



AG(V) = -p.\Afi\V + 6a pp V 2/3 + 12v pp V 1/3 + k B T\n ^6 (J^J ' \ + const. (4.6) 



with <j pp = ctioo- Similarly to the isotropic case, in the cluster free energy (|4.6|) both a 



logarithmic term and an offset are added to the classical CNT expression of AG for a cubic 



cluster of side V 1 ^ 3 . The Tolman term in Eq. (4.6) only appears if we envisage an energy 
penalty, that is v pp per unit length, also for the edges. 

More generally, in all the anisotropic-nucleation models examined in Appendix D, the 
consideration of clusters of same type but unequal edges/semiaxes provides for "breathing" 
fluctuations of the surface that determine the appearance of a logarithmic term in AG. In 
fact, for all such models, the analytically computed AG(V^) is asymptotically given, as in 



Eq. (4.6), by the CNT expression - as written for the respective symmetric shape - plus 
subleading terms in the form of a Tolman term, a universal logarithm (c/c^T In V^ d ~ 1 ^ d in 
d dimensions), and a negative offset. The value of c is 1/2 for rectangles and 1 for both 
cuboids and ellipsoids. This is to be contrasted with the quasispherical-cluster case, where 



c = —7/6 by Eq. (4.4 ) . Apparently, the value of c is sensitive to both the space dimensionality 
and the number of independent parameters that are needed to describe the cluster shape, 
in turn crucial to determine the entropy contents of the surface degrees of freedom (for a 
quasispherical cluster, this number of parameters goes to infinity with V). In short, a large 
anisotropy in the interface tension has the overall effect of drastically reducing the spectrum 
of thermal fluctuations of cluster shape. The reduction cancels the entropy gain which these 
fluctuations produced in the isotropic case. 

This attractive prediction is a difficult one to fully validate numerically at present. A 
logarithmic correction to CNT can only be detected if we push the numerical investigation 
of AG(V) so close to coexistence as to make the Dillmann-Meier form exact for all but the 
smallest clusters, and that is still a difficult task (see more in the next Section). In the near 
future, with faster computers becoming available, we can imagine that it will be possible to 
directly probe the cluster geometry through the optimization of the logarithmic prefactor 



in an ansatz of the kind (4.4) or (4.6), and thus choose among the many cluster models on 



the market the one which is most appropriate to the problem at hand. 
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V. NUMERICAL ASSESSMENT OF THE THEORY 



We now critically consider if there are signatures of the degree of anisotropy of the 
interface free energy in the free-energy cost of cluster formation for a specific instance of 
microscopic interaction. 

We first recall how the work of formation of a n-particle cluster is calculated from sim- 
ulations! J32H3I]- Given a criterion to identify solid-like clusters within a predominantly 
liquid system of N particles, the average number of n-clusters is given, for 1 <C n ~ n*, by 
N n = N e - ^ Gn-nw ), where \i\ is the chemical potential of the liquid and G n is the 0(n) Gibbs 
free energy of the n-cluster, including also the contribution associated with the wiggling of 
the cluster center of mass within a cavity of volume V/N (observe that CNT estimates G n as 
njj, s + ccr(n/p s ) 2 / 3 , where jx s is the chemical potential of the solid and c a geometrical factor). 
For rare clusters, it thus follows that AG(n) = G n — n/ii = —kBT\n(N n /N). This equation 
is then taken to represent the work of cluster formation for all n > 1. Maibaum [35J has 
shown that the same formula applies for the Ising model. 

However, for quenches that are not too deep, the spontaneous occurrence of a large solid 
cluster in the metastable liquid is a rare event. This poses a problem of poor statistics 
in the Monte Carlo (MC) estimation of N n , which is overcome through e.g. the use of a 
biasing potential that couples with the size n max of the largest cluster. In practice, this keeps 
the system in the metastable state for all the n's of interest. By properly reweighting the 
sampled microstates one eventually recovers the ordinary ensemble averages. This umbrella- 
sampling (US) method was used in Refs. [32} [36] to compute AG(n) for the Lennard- Jones 
fluid and the 3D Ising model, respectively. The main obstacle to the calculation of AG(n) 
by US is the necessity of performing the identification of the largest cluster in the system 
after every MC move. This problem can be somewhat mitigated by the use of a hybrid MC 
algorithm [37], which in our case reduced the simulation time by a factor of about 20. 

A low-temperature Ising magnet where the majority of spins point against the applied 
field probably yields the simplest possible setup for the study of nucleation. Along the 
first-order transition line of the model, where two ( "up" and "down" ) ferromagnetic phases 
coexist, the interface (say, (100)) between the two phases undergoes a roughening transition 
at a certain T = Tr. The up-down interface tension at coexistence is strongly anisotropic 
close to zero temperature; moreover, it is either singular or smooth according to whether T 
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is below or above Tr. Strictly speaking, the interface tension is anisotropic also above Tr, 
though less and less so when approaching the critical temperature T c from below [38, 39J. 
Exactly at T c the interface tension critically vanishes (10] • When a sample originally prepared 
in the "down" phase is slightly pushed away from coexistence by a small positive field and 
thus made metastable, the critical droplet of the "up" phase is expected to be less and less 
spherical as T decreases. 

With the 3D Ising model as a test system, we carried out a series of extensive US sim- 
ulations, computing the cluster free energy AG(n) relative to the nucleation process of 
magnetization reversal for a fixed T = 0.6 T c , slightly above the roughening temperature Tr 
of the (100) facet (Tr = 0.5438 . . . T c [S]), and for a number of values of the external field h 
(0.30, 0.35, . . . , 0.65, in J units). Two up spins are said to belong to the same cluster if there 
is a sequence of neighboring up spins between them; the counting of clusters was done with 
the Hoshen-Kopelman algorithm jl2]. The absolute value of AG(1) was determined through 
a standard MC simulation of the system with all spins down, with no bias imposed on the 
sampling of the equilibrium distribution. We point out that, at the chosen temperature, 
the Ising surface tension is barely anisotropic |38j , which would exclude a net preference for 
either the spherical or the cubic shape. Furthermore, we are sufficiently far away from T c 
not to worry about the percolation transition of geometric clusters which was first described 
in jl3]. This event, which would invalidate the assumption (at the heart of the conventional 
picture of nucleation) of a dilute gas of clusters, is still far away here. 

Coherently with the physical picture at the basis of our theory, we verified for all the h 
considered that clusters close to critical indeed contain the vast majority of up spins in the 
system. A sample of the critical cluster for h = 0.30 is shown in Fig. 1. Looking at this 
picture, it is hard to say whether this particular realization of the critical cluster resembles 
more a sphere or a cube. When moving to h = 0, a spherical shape is eventually preferred 
over the cube far above Tr, whereas the opposite occurs much below Tr. 

In Fig. 2, the ratio a l {n) of the surface free energy F s {n) = AG(n) + \ Afi\n to the area 
S(n) of the cluster surface is reported as a function of n" 1 / 3 , and the data are fit using 



the functions (4.4) and (4.6) (we stress that different expressions apply for S(n) on the left 
and right panels of Fig. 2, i.e., (367r) 1 / 3 n 2 / 3 and 6n 2//3 respectively; accordingly, the spherical 
cr's would typically turn out a factor 6/(367r) 1 / 3 larger than the cubic cr's). Both fits are 



based on three parameters, namely a, 5, and e, which enter in a different way in Eqs. (4.4) 
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FIG. 1: (Color online). A snapshot taken from our Monte Carlo simulation of the 3D Ising model 
at T = 0.6 T c and h = 0.30, showing a cluster of n = 685 up spins, i.e., close to the critical size for 
that h. Up spins are differently colored according to the number of nearest-neighboring up spins 
(blue, 6; cyan, 5; green, 2-4; magenta, 1; red, 0). Down spins are not shown. 



and (4.6). However, the dependence on n is similar for the two fitting functions, except for 
the numerical factor in front of the (parameter-free) n~ 2 ' 3 lnn term. Looking at Fig. 2, it 
appears that the quality of the "cubic" fit is slightly better than that of the "spherical" fit, 
in line with the fact that, for T > Tr, the Ising surface tension is moderately anisotropic. 
Clearly, at T = 0.6 T c the nucleus is neither spherical nor cubic, and one may object that 
neither of the fits would actually be meaningful. We nonetheless argue that, within the 
uncertainty associated with the finite h value in the simulations, the better one of the fits 
will correspond to the regular shape which is closest to that of the real nucleus, thus giving 
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FIG. 2: (Color online). The cluster free energy a 1 of the 3D Ising model on a cubic lattice in 
units of J /a? is plotted as a function of n -1 / 3 (and up to 80 -1 / 3 ) for three values of h and for 
T = 0.6 T c (a is the lattice spacing and J > is the spin-coupling constant). The lattice includes 
20 3 sites (25 3 for h = 0.30). Umbrella-sampling simulations consisted of 4M equilibrium sweeps 
for each n window (one window covering eleven values of n). Thick colored lines, MC data; black 



lines, least-square fits of the n > 80 data points for h = 0.30,0.40,0.50, based on Eq. (4.4) (left) 



and (4.6) (right). Data plotted in the two panels look different simply because the expressions of 
cluster area S(n) are different between left an ight (see text). Inset, the difference between the 
raw data and the fit. 



a qualitative indication of the prevailing isotropic or anisotropic character of the solid-liquid 
interface tension. When going to smaller and smaller h, and provided T is sufficiently above 
T R , we expect that the "spherical" fit would eventually become better than the "cubic" fit. 

VI. CONCLUSIONS 

In order to estimate from nucleation the solid-liquid interface free energy a m of a sub- 
stance, two indirect routes are available: one is through the measurement of the solid nu- 
cleation rate as a function of temperature (see e.g. [161 - H8] ). the other is via the free energy 
AG of solid-cluster formation in a supercooled-liquid host, as determined for example in a 
numerical simulation experiment for a system model. In both cases, the theoretical frame- 
work of classical nucleation theory (CNT) has routinely been employed to extract a m . This 
is far from satisfactory, as discussed at length in Ref. [Sj and in many other papers, due to 
the neglected cluster interface-tension dependence on both the droplet volume V and the 
supersaturation |A/i|. 

Concentrating on the expression of the cluster formation energy AG as a function of 
V and Afi, we gave here an extension of the modified CNT theory first introduced in [8], 
now including anisotropy, which is important when only a few interface orientations survive 
in the equilibrium average cluster shape. We showed that, also in this universal 
non-CNT In V term is found in the asymptotic expression of the surface free energy versus 
volume, so long as an infinity of regular shapes is allowed to occur. However that term 
has now a different prefactor with respect to the quasispherical case. In particular, the 
sign is positive for large anisotropy and negative for vanishing anisotropy. The sign of 
that prefactor, which we surmise is related to the amount of surface entropy developed by 
cluster shape fluctuations, is proposed as the imprinted signature of the geometrical shapes 
most preferred by the nucleation cluster - negative for spherical or very isotropic shapes, 
positive for nearly polyhedric or anyway very anisotropic shapes. For the 3D Ising model 
slightly above the (100) roughening transition temperature, the detected sign suggests cubic 
rather than spherical cluster symmetry for moderate supersaturation/external field. Much 
more work and larger simulation sizes should be needed in the future in order to verify the 
expected change of sign of the In V term as spherical shapes will be approached closer and 
closer to the coexistence line when the temperature is quite larger than Tr (though still far 
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from the critical region). 
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Appendix A: Calculation of a m and 5 n 



In this Appendix, we provide approximate expressions for the quantities a m and 5 m in 



Eqs. (2.8) for a specific model of homogeneous-system free energy g(<p) in the functional 



(2.2) 



Once the exact OP profile (f>o(z) of the planar interface has been determined for the given 
g, the explicit values of a m and S m , and of e m = e L (T m ) can be computed. While a m and 
e m are strictly positive quantities, the sign of S m is not a priori definite. A special but 
sufficiently general case of g function is the following: 



g{<j>; T = T n 



C20<P 



1 + (75 + 27 6 ) — + 7e- 



4 



(A.l) 



with 75 > — 1 — 37 6 for < 7 6 < 1 and 75 > — 27 6 — 2^/Te" for 7 6 > 1. Equation (A.l) is 



the most general sixth-degree polynomial which admits two non-equivalent minimum valleys 
at and s o, and no further negative minimum between them. For this g, the differential 
equation (2.4) is still too difficult to solve in closed form for generic k, even when 75 = 76 



0. 



Hence, we decided to work perturbatively in k, 75, and 76. 

At zeroth order, i.e., k = 7.5 = 76 = 0, corresponding to 4 theory, the solution to (2.4) is 

'z-C s 



75 = 76 

M z ) = 



■PsO 

2 



tanh 



(A.2) 



with £ = a/2c/c2o and arbitrary C. We fix C by requiring that the interface is centered at 
z = (hence C = 0). Then, by still keeping 75 = 76 = 0, we switch on k and search for a 



second-order solution to Eq. (2.4) in the form 



<Pq(z) = 0oO) + ^pXi(z) + X2(z). 



(A.3) 
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/ — "/ 1 —112 



We thus arrive at the two equations: 

c&'aX'i - 9o&o)Xi = cf ( 0o0o - 7$o 

and 



where 



9o{4>) = C2O0 2 1 



By requiring that <fio(z) is centered at z = we obtain 



>s0 



and 

X2{z) 



>s0 



cosh 2 (z/£) 



32 tanh d 



cosh 2 (^/£) 



12- tanh 2 



2tanh 



z z 



(A.4) 



(A.5) 



(A.6) 



(A.7) 



; tanh - + 2 



tanh 



z _z 



(A.8) 



Hence, we find 5 m = since the function c0' o 2 (z) + 2k0 / o /2 (z) is even. Actually, the result 
5 m = is valid at any order in k when 75 = 7 6 = (see below). Up to second order in k, 
the values of a m and e m are given by: 



2 k 38 / k \ 2 
1 + 5 ^2 ~ 35 Vd^J 



c0; 



s0 . 



3£ 



tt 2 - 6 /26 7T 2 \ « /1566 4tt 2 \ / k \ 2 



12 



5 



c£ 2 + V 



175 



(—)' 
\cPJ 



(A.9) 



Next, we take k, 75, and 76 all non-zero and of the same order of magnitude, and search 



for a first-order solution to (2.4) in the form 



M z ) = 0oO) + 75^1(2) + 76&(z) + ^2^i( z ) • 



(A.10) 



Upon inserting (A. 10) into Eq. (2.4), we obtain two independent equations for ipi(z) and 
£1(2), namely 

cM-g'oM^-^^ 1 (A.n) 

and 



(A.12) 



^0 </^o 
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while Xi{ z ) is still given by Eq. (A.7). The solutions to Eqs. (A.ll) and (A. 12) such that each 
term of (A. 10) separately meets the requirement of being centered at zero are the following: 



',cosh\z/£) 



^1 — In 2 + - — In cosh 



and 



^•0 



) cosh 2 (z/i) 



z z I z 
3(1 — In 2) + 3- — 3 In cosh tanh - 

•C K, 2 -C 



(A.13) 



(A.14) 



Upon plugging the by now specified (fio(z) in the integrals defining a m ,5 m , and e m , we 
eventually obtain the formulae: 



, 1 13 2 « \ ccj) 2 s0 
am = (l + 4 75 + 20 76 + 5 ^JlF' "'" 



S, 



48 



7T 2 — 6 7T 2 — 6 



12 



48 



"75 



17vr 2 1 



240 



76 



(75 + 37e ) i , and 

26 7T 2 \ K 



(A.15) 



We thus see that 5 m is generically non-zero and may be of both signs. 
In conclusion, we give a proof that 5 m vanishes identically for 

V 4>so 



(A.16) 



whatever k is (a different argument can be found in [13])- Let </>(z) be a solution to Eq. (2.4) 
obeying the boundary conditions 



-00 



'+00) = 0, d/(±oo) = <b"(±oo) 



0. 



(A.17) 



There is an infinite number of such solutions, differing from each other by a simple transla- 
tion. Let us first prove that <j)(z) = (f) s o — 



-z) is also a solution to (2.4). We have: 



g(<K*)) = g(<K-z)) ; </>'(*) = </>'(-*); = -<?{-*); A*) = A-*) ■ (A.is) 



We thus see that 



K0{-z)r{-z) - C ^\-z) - ^" 2 {-z) + g{<t>{-z)) = 



(A.19) 



since Eq. (2.4) is satisfied by <fi for any z. Hence, <fi(z) obeys the differential equation (2.4) 



Moreover, like <p(z), <p(z) also satisfies the conditions (A.17). This is not yet sufficient to 



conclude that <p(z) and <fi(z) are the same function since they could differ by a translation 
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along z. However, if among the infinite possibilities the one is selected such that 0(0) = 
</> s o/2, then 0(0) = 4> s q/2 and the two functions coincide: <f)(z) = <fi(z), implying 



z)+<b-z) 



>s0 



for any z . 



(A.20) 



Upon differentiating (A.20) with respect to z we find that <p'(—z) = 4>'(z) and the function 
•(z) is even. This is enough to conclude that j dz z<fi'{z) = (the interface is centered in 0). 

—z) = —<f)"(z) and 4>"{z) is an odd function 

of z (while 



Differentiating (A.20) once more, we obtain </ 
l )" 2 {z) is even). As a result, a m 5. 



- / dzz[c(p ,2 (z) + 2k<P" 2 (z)} = and the 



proof is complete. 



Appendix B: Derivation of Eq. (3.4) 



Let i be a curve in E parametrized by the arc length s and denote (t, n, b) the Frenet 
trihedron in H(uq,Vq) € £. Note that we are using nearly the same symbol for the normal to 
E (n) and for the normal vector to I in R(mo,i>o) (n), though the two vectors are generally 
distinct. Now consider the Darboux frame (T,N,B) with T = t, N = n (the unit normal 
to E in R(uo, ^o))> an d B = T A N. Clearly, by a convenient rotation around T = t, n and 
b are carried to N and B, respectively. Calling a(s) the rotation angle, 



N 



/ 



V 



i o o 

cos a sin a 
— sin a cos a 



\ 



n 



Using the Frenet-Serret formulae, namely 

dt 

d~s ~ 
dn 

d7 = 

db 

d7 " 



/tn : 



-K,t 



rn 



(B.l) 



(B.2) 



where k is the curvature and r is the torsion of £, we easily get 



' dT/ds^ 

dN/ds 
\^dB/ds J 



( 







\ 



-Kg -Tg 



N 



(B.3) 
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where n n = ncosa is the normal curvature, n g = — /tsina the geodetic curvature, and 
r g = t + da/ds the geodetic torsion. 

For the it-lines, if we identify T with R u then B = T A N = R, Similarly, for the 
v -lines, if we identify T with R v then B = T A N = R„. We thus obtain: 



dn 

du 



'^n ^-f^-w i~g Rv and — — Tg ^R u ^R. v . 



(B.4) 



Moreover, 



dR u 

du 
dR v 

du 



k,^R v + K^h and 



dRu 

dv 



~ K g Rv T g ^6 I 



dR,, 



+ and = k^R u + K^h . 



(B.5) 



From R,, 



R OT( , we derive 



4 1 )R M + 4 2 )R,+ (r« + rf)n = 0. 



(B.6) 



Since R u , R v , and n are linearly independent, it necessarily follows that 



K W = K £) = Q and T m = _ r W 



- T 9- 



(B.7) 



Since the geodetic curvature vanishes, our assumption that R u and R„ are orthonormal 
vectors implies that the coordinate lines are surface geodesies. 



Finally, putting Eqs. (B.4) and (B.7) together we promptly get Eq. (3.4). 



Appendix C: Small fluctuations about a spherical interface 



We here provide the detailed derivation of Eq. (4.4) for the free energy of a quasispherical 



interface E. The starting point is the expansion of the relative amount of asphericity, e( 



in real spherical harmonics, Eq. (4.3). In view of the smallness of the expansion coefficients 
Xi >rn , the enclosed volume and area of S can be approximated as 



V[E] 



dSr ■ n = ^Rl + R 3 J2 x lm = /({*}) 

l>0,m 



and 



A[E] = j dS = AttRI + £ p + i + 2 ) x\ m ee 4kB* g({x}) 

^ S l>0,m 



(C.l) 



(C.2) 
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f({x}) and g({x}) being close-to-1 factors. In writing the two formulae above we supposed 
x ,o = 0, which can always be assumed by suitably redefining in R(9, </>) the radius Rq and 
the other coefficients x^ m . In order to evaluate the mean curvature H, we start from 

V .. =J _£(^ i a i 

r A or r sin 6 oa r sin U o<p 



where 



i 1 2 1 C ' 

n r = 1 — -e fl — 



2~° 2 S m 2 e , 
he = -eg(l - e) ; 

ft, = -^^. (C.4) 
smv 



From that we get 



V ' fi = R(hj i 1 + \ L2<6 ' 0) ~ 2^' 0) ) ' (C - 5) 

where 

1 9 / . J\ 1 9 



sinfldfl V 90/ sin 2 
Eventually, we obtain: 

I dS (a m - 2a m 5 m H + 2XH 2 ) = Ana m R 2 + ^ (/ 2 + / + 2) 



X 



Lm 



-8ira m 5 m R - a m 5 m R K l + l ) x lm + §7rA + \ Yl l ( l + X )( Z " X )( Z + 2 Kr. 

£>0,m Z>l,m 



(C.7) 



Finally, we specify the integral measure in (4.2): 



^=/ + ~n (!<o r^, (c.8) 



where 5 = (36vr) 1 / 3 1/ 2 / 3 is the area of the spherical surface of volume V and s = Ana . 



Equation (C.8) follows from requiring that the present theory (in fact the theory with an 
upper cutoff on I, see below) should coincide with the continuum limit of the field theory 
for a solid-on-solid (SOS) model with real heights defined on nodes uniformly placed over a 
sphere of radius a/ ~S] (4~7r). 
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To prove this, first observe that the equation for the generic £ entering in the functional 



integral is R - R = Yli>o, m R Y hm {9, <f>)xi >m . Since R = ^S/(Att) up to terms 0(x 2 m ), the 
height profile which the equation for E corresponds to is 

~S 



l>0.m 




4tt 



Yi, m (VLi)xi. 



(C.9) 



for i — 1, 2, . . . , n and n = (l max + l) 2 — 1 — S/a 2 (the necessity of an upper cutoff l r 



on 



/ given by the following Eq. (C.17) will be motivated later). The relation between the two 
theories passes through the identification 



n 

i=l 



dhi 



< — > 



J 



n ^ xi ' 



(CIO) 



l>0,m 



where J is the Jacobian of the transformation (|C.9|): 

dh„ 



J = det 



dx 



I,m 



g \ n/2 
47T 



max max 



(C.ll) 



Yi-i(Q„) ... ^ max ,/ max (^n) 

Called AQ = Air/n the element of solid angle assigned to each node we have 



Yi :m {p,i)Yii >m i {Q,i 



AQ 



d 2 nY ljm (n)Y rtm/ {fi) 



AQ 



(C.12) 



Hence, for sufficiently large n the columns of the matrix (C.ll) are mutually orthogonal 



n-vectors. In order that every column vector be normalized, it suffices to multiply the 



whole matrix in (C.ll) by \/Af2, thus getting an orthogonal matrix (of unit determinant). 
Therefore, we find 

which amounts to take s = Ana 2 in Eq. (C.8). This completes our proof. 



We can now go on to compute the partition function (4.2). We first calculate the integral 
on R by rearranging the delta function in Z s as 



S[-nRlf({x})-V 



s(R -[Anf({x})/(3V)]- 1/3 



(C.14) 



(367r) 1 /3F2/3J({ x })l/3 

After doing the trivial integral on Rq, we remain with a factor /({x}) -1 / 3 which, within a 
quadratic theory, can be treated as follows: 

-1/3 



/(to)' 



-1/3 



L + 4?r Y X lm 

l>0,m 



47r Y X l™ - eXP ] 4vr Y X l™ f 
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(C.15) 
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In the end, we arrive at a Gaussian integral which is readily computed: 



(36tt) 



-1/3 
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exp < /3p s \Afi\V — /3a m S — 87T/3A + %Kf3a m 8 rt 



3V 

47T 



1/3' 



Z>1 



S n- 1 

2^S. 



(3cr m S , 



1 + l ^f-(l 2 + l-2) + 2tt/3XI(1 + - 1)(Z + 2) 
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(C.16) 



Without a proper ultraviolet cutoff Z max the / sum in lnZ s does not converge. This is a 
typical occurrence for field theories on the continuum, which do not consider the granularity 
of matter at the most fundamental level. We fix Z max by requiring that the total number of 
(l,m) modes be equal to the average number of SOS heights/atoms on the cluster surface. 
It thus follows: 



S 



1 . 



(C.17) 



With this cutoff, the surface free energy becomes F s = a(S)S, with an interface tension 
cr(S) dressed by thermal fluctuations: 

y/S/a-1 



a (S)=a m + k ^- (2l + l)\n[A + B(l 2 + l-2) + C(l 2 + l-2) 2 } 

1=2 
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4tt 
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(C.18) 



The quantities A,B, and C in Eq. (C.18) are given by 



A = — R = 2C ° i D ° i — r= — 
S 2 ' S 2 + sVs S ' S 2 ' 



where 
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(C.20) 



By the Euler-Mac Laurin formula, the residual sum in Eq. (C.18 ) can be evaluated explicitly. 
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After a tedious and rather lengthy derivation, we obtain (for A 7^ 0): 
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(C.21) 



up to terms o(S f : ). We wrote a computer code to evaluate the sum in (C.18) numerically 



for large S, and so checked that every single term in Eq. (C.21) is indeed correct. 



Appendix D: Anisotropic-interface models of nucleation 

We here show that nonperturbative corrections to CNT do also arise when the interface 
tension is infinitely anisotropic. In this case, the admissible cluster shapes are all regular 



and the functional integral (4.2) is greatly simplified, reducing to a standard integral over 



the few independent variables which concur to define the allowed clusters. The terms in 



(3.15) beyond the surface-tension term do also contribute to the total surface free energy if 



anisotropy is strong but not infinitely so. 



Our argument goes as follows. Let us, for instance, consider the interface tension (4.5). 



For M 1, we expect that the leading contribution to the functional integral (4.2) be 
given by rectangular cuboids with slightly rounded edges and vertices. Since we are only 
interested in making a rough estimation of the relative magnitude of each contribution to 
H s , we assume that the surface of a rounded edge is one fourth of a cylindrical surface 
(H = 1/a) whereas that of a rounded vertex is an octant of a sphere (H = 2 /a), a being 
a microscopic diameter. Also observe that: aioo ~ fcsT m /a 2 ; the average value of er(n) on 
an edge or vertex is ~ Mciioo; the Tolman length is 5 m ~ a; A is roughly k/c times a m , 



hence A ~ a m a . We now decompose (3.15) into the sum of three integrals, respectively over 



faces, edges, and vertices. Denoting li,l2,h (all much larger than a) the side lengths of the 
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cuboid if its edges and vertices were taken to be sharp, the integral over faces is practically 
equal to 2<7 10 o(M2 + his + hh) ~ &B^m(V2 + hh + hh)/® 2 ] U P to a factor of order one, the 
integral over edges is given by M/c j eT m (/ 1 + l 2 + h)/ a ] finally, the integral over vertices is of 
the order of MfceT m . We then see that, when M — > 00 for fixed a, only faces contribute to 



the integral (4.2), while edges and vertices would only matter if M were finite. 

In general terms, from the knowledge of the "Wulff plot" cr(n), the equilibrium cluster 
shape follows from the so-called Wulff construction |44j : (i) draw the planes perpendicular 
to the unit vectors n and at a distance er(n) away from the origin; (ii) for each plane, discard 
the half-space of M. 3 that lies on the far side of the plane from the origin. The convex region 
consisting of the intersection of the retained half-spaces is the cluster of lowest surface energy. 
When cr(n) is smooth, this Wulff cluster is bounded by part of the envelope of the planes; 
the parts of the envelope not bounding the convex body - the "ears" or "swallowtails" which 
are e.g. visible in Figs. 3 and 6 below - are unphysical. 



1. Rectangles 

In two dimensions, a rudimentary model of nucleation is that which only allows for 
rectangular clusters. This is relevant for two-dimensional crystals of square symmetry, and 
could be obtained from a smooth interface tension of the form 



a 



(0) =<7io [l + Msin 2 (20)] (D.l) 



upon taking the infinite-M limit. In Eq. ( |D.l ), <p is the polar angle of the normal vector 



while 0io is the free energy of the cheapest, (10) facet. As M — > 00, all normal directions 
different from [10], [01], [10], and [01] are excluded from the equilibrium cluster shape and a 
perfectly square surface is obtained. This is illustrated in Fig. 3 for three values of M; here 
and elsewhere, the envelope of perpendicular planes is given, in parametric terms, by the 
equations [45] : 

x = cos (j) cr(0) — sin <pa'((f)) ; 

y = sin 0cr((/>) + cos (f) a' ((f)) . (D.2) 



Note that the edge fluctuations deforming the square in a rectangle are still allowed by 
(D.l ) in the infinite-M limit, since rectangles and cubes share the same type of facets. Hence, 
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FIG. 3: (Color online). Two-dimensional Wulff construction for a square cluster. Left: Polar plot 



of the interface tension a in Eq. (D.l) (red curves) as a function of the normal n = cos0x + sin<^y 
to a cluster face, for three distinct values of M. Right: Equilibrium cluster shape as the envelope 
of the family of perpendicular planes (blue curves). In this case, the cluster of minimum surface 
energy (colored in cyan) has curved faces, but sharp corners. The envelope continues beyond the 
corners, but these parts have no physical meaning. 

assuming that only rectangular shapes have a non-zero Boltzmann weight in the functional 



integral (4.2), the surface free energy reduces to: 

"+OO 



-i r-\-co p-voo 

F s (V) = --\n da dbe~ 2l5a{a+b) 5(ab-V) (D.3) 
P Jo Jo 
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FIG. 4: (Color online). Rectangles: (3AG(V) vs. V for fig = j3(X = 1 (black solid line). The dotted 



red line corresponds to the approximant (D.8). Inset, the difference between (D.8) and /3AG(V) 



We easily find: 



0AG(V) = -(3gV + f3F s {V) 

/+oo 
dx exp {-r(coshx - 1)} 
-oo 



-PgV + AfiaVV - In \ 2e 



+oo 



-Tt 



dt 



-PgV + 4f3aW - In {2e T K {r)} , 



(DA) 

where g is proportional to the supersaturation, r = 4f3a\/V, and Kq is a modified Bessel 
function of the second kind. The last term in Eq. ( |D.4 ) is the full correction to CNT as 
formulated for squares. A typical profile of (3AG(V) is plotted in Fig. 4. 

At variance with CNT, AG(V) shows a weak divergence to — oo for V — > 0, due to the 
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absence of a lower cutoff volume. For r < 1, 



K q (t) = -ln(r/2) - 7 + 0(r 2 In r) , (D.5) 

with 7 = 0.5772 . . . (Eulero-Mascheroni constant). Hence, the singular behavior of /3AG(V) 
for small V is of the kind 



(3AG{V) ~ - In | - In feaVv) } . (D.6) 
Conversely, for large r values, 




A/^ (D.7) 



and we obtain 



PAG(V) ~ -/fyV + 4/9(tVF + i In - - hi(2vr) . (D.8) 



The goodness of the approximation (D.8 ) can be judged from the inset of Fig. 4, which shows 



that the approximation is accurate for all values of V but for the smallest ones. 
2. Truncated rectangles 

In order to study the effects on nucleation of a more complicate type of interface-tension 
anisotropy, we further enrich our book of patterns, passing from rectangles to truncated 
rectangles. By the name of truncated rectangle we mean the octagon represented in Fig. 5. 
This occurs when the cost of (11) and equivalent facets is of the same order of crio, while all 
other facets are much higher in energy and can be ruled out. 

A Wulff plot giving origin to truncated squares is: 



<r(4>) = °"io 



1 sin 2 (20) + M sin 2 (40) 

0"io 



(D.9) 



with infinite M (see Fig. 6). The polar plot of (D.9) for finite M is a smoothed eight-pointed 



star with hollows at the normal directions satisfying sin(40) = 0. Depending on the ratio 
of (Tii to a w , the equilibrium cluster shape shows (i) just (11) facets (crn/cio < v2/2); (ii) 
both (11) and (10) facets (y/2/2 < a u /a 10 < ^2); (iii) just (10) facets (on/trio > V2). 

In order to prove this, we observe that, for fixed a, b, and £ (with £ < £ max = 
(l/\/2) min{a, b}), the "volume" and "area" of the truncated rectangle are given respec- 
tively by V = ab - £ 2 and A = A u + A w , with A n = A£ and A 10 = 2(a + b- 2y/2£), leading 

30 









V 


\ 




h 5 a H 



FIG. 5: A truncated rectangle. For fixed a and b, the maximum £ value is £ max = (l/\/2) min{a, b}. 
to a surface energy of 

E s = Aa lx l + 2cr w (a + b - 2y/2£) . (D.10) 

To determine the cluster shape at zero temperature, E s should be minimized as a function 
of a, b, and £ under the constraint of a fixed ab — i 2 (= V). Setting a = xy/V and b = y\/V 
(with x, y > 0), we are led to minimize 4(<7n — \/2<Ji )\/xy — 1 + 2<7 10 (:r + y) as a function 
of x and y. By a straightforward calculation we find: 

x = y=l {£ = 0), for — > V2; 

ClO 

1 / * V2-a u /a 10 \ V2 a u 

x = y = = — = = = , for — < — < V 2 : 



V / (710 2 



s = y = V2 (^L = lV for^<^. (D.ll) 
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FIG. 6: (Color online). Two-dimensional Wulff construction for the a model at Eq. (D.9), with 
M = 1. Left: Wulff plot (red) for three values of on/cxio- Right: Equilibrium cluster shape (the 
boundary of the cyan-colored region). The parts of the envelope of the family of perpendicular 
planes beyond the corners of the cluster are unphysical. 

The equilibrium cluster shape is then a square (respectively, a 45-degree tilted square) for 
°"n/°"io values larger than \/2 (smaller than y/2/2), while being a truncated square otherwise 
(see Fig. 6). 

Now going to the nucleation model for truncated rectangles, the surface free energy reads: 

F S (V) = — In / / dadb / ^ e -4^ 11 ^-2^ 10 (a + 6-2^) 6 , ab -f-v), (D.12) 
P J Jo Jo to 
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where £o is an arbitrary length. By integrating the delta out, we obtain: 
/3AG(V) = -f3gV + 4f3a 10 VV-Un(^j 

/+oo r+oo ^_ 

dx dyQ (min{e x ', (1 + y 2 )e" x } - V2y 

x exp{-(r 11 -v / 2r 10 ) 2 /}exp{-^(e a + (l+y 2 )e^-2)} , (D.13) 
where is Heaviside's function, t w = A/3aioVV, and m = 4/3<JuVV. For (3g\q = 1,£$ = 1) 



and <j nl 0io = 0.5, 1, 2, 20, 200, the plot of (D.13) is reported in Fig. 7. Note that a precritical 



minimum shows up for any finite value of 011/010, which moves toward zero upon increasing 
the interface-tension anisotropy. An even more complex behavior is seen for 0h/0iq = 0.5, 
where a bump emerges beyond the critical maximum. 

When an ^ 0io, it is natural to expect that the model of truncated rectangles reduces 



to the rectangular-cluster model. This can be proved analytically starting from Eq. (D.12). 
First, a, b, and £ are rescaled by dividing by \/V ; then one observes that 

^/y^anWt „ -J—5(t) . (D.14) 

2/3011 

Hence, aside from a constant equal to ln(4/?0n), the f3AG(V) function for truncated rect- 
angles merges, for very large 011/010, into the analogous function for rectangles. This fact 
is shown numerically in the inset of Fig. 7. 

3. Rectangular cuboids 



When the Wulff plot is as in Eq. (4.5) with infinite M, the only admissible shapes are 
rectangular cuboids. Denoting a, b, and c the edges of a cuboid, the V^-dependent surface 
free energy is defined as 

[• r p+co 

pF s (V) = -In jjj dadbdce' 2 ^ iab+ac+bc) 5{abc-V) 

= -ln/^ + ~d a d6ie X p{-^( 6 + °± 6 )}, (D.15) 
where r = 6(3aV 2 ^ 3 . With another change of variables, we arrive at 



PAG(V) = -PgV + QfiaV 2/3 -\n j! + dx dy exp { - - (e x+y + e~ x + e~ y - 3) } . (D.16) 



33 




FIG. 7: (Color online). Truncated rectangles: (3AG(V) vs. V for f3g = (3a\o = 1 e 4 = 1. A 
constant of ln(4/3<7n) was subtracted from (3AG(V) in order to garantee the confluence of its plot 
to that for rectangles, in the limit crn/cio — > +oo. A number of an/cio values are considered: 
0.5 (black), 1 (blue), 2 (cyan), 20 (magenta), and 200 (red, practically indistinguishable from the 
rectangular case). In the inset, we zoom on the small- F region, evidencing the singular behavior 
of AG(V) for V — > 0. Apparently, for all finite an values, the curve blows up to +oo rather than 
to — oo, as instead occurs for rectangles. 

The above formula is well suited for the numerical evaluation of AG(V). For f3g = /3a = 1, 
the profile of /3AG(V) is plotted in Fig. 8. 

In order to discover the analytic behavior of AG(V) at small and at large V's, we should 



further elaborate on Eq. (D.16). Setting a + b = x and ab = y in (D.15), a and b are the 



solutions to the equation t 2 — xt + y = 0, whose discriminant is non- negative for x > 2y/y. 
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FIG. 8: (Color online). Rectangular cuboids: /3AG(V) vs. V for (5g = f3a = 1 (black solid line). 
In blue, the same function when we include the cost of the edges {f3g = (3a = f3v = 1). The red 
dotted line is the approximant (D.27) while the cyan dotted line is the approximant ( |D.31 ). In the 



inset, we plot in red the difference between the approximation (D.27) and (D.16), and in blue the 



difference between (D.31) and (D.28). 



Moreover, the Jacobian of the transformation is 1/yi 2 - %. Hence, we get f3F s (V) 
-\nI{V) with 

(D.17) 



1=2 r°dy r^M-mk±m 

Jo J2^y y^/x 2 -4y 



With the further transformations x — > z = \/ x 2 — Ay and 



z — > w 



z/y, we eventually 
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obtain: 



o ' y Jo ' v^ 2 + A v 

-+oo e -(r/3)y r+oo e xp|-Ww 2 +^ 

2 1 dy / dw 

h y Jo ■ / ™ 2 + i£ 



Since 



we finally find: 



+ °° , expi-Vx 2 + c 2 } 
dx 



o 



-oo _i 
dt 



2 _ ~2 



^o(c) 



1 = 2/ dy 6 - K (f- 

Jo y \^Vy 



j 4t 3 1 tf (x) 
4/ d ^ ex P^-^2-j 



(D.18) 



(D.19) 



(D.20) 



In Eq. (D.20) we recognize a particular Meijer function, Gq®((t/3) 3 \0, 0, 0), whose behavior 
at small r is: 
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(In r) 2 + 9(7 - In 3) In r + 0(1) 



(D.21) 



From the above, we can draw the main singular term in /3AG(V) at small V, that is 



/3AG(V) ~ -2 In (-ln(6/3ay 2/3 )) 



(D.22) 



which is similar to (D.6). 



The large- y behavior of AG(V) can also be obtained from Eq. (D.20). For r > 1, we 



are allowed to replace K (x) with Eq. (D.7) and thus estimate I through the integral 



-+oc exp <j -x - 
= 2^2^ / dx ., 2 

o x 



[D.23) 



Suspecting a dominant term of r in — In 1^, we consider 
e T 'I OB = 2y^-J d2z- 3/2 exp|r (\ 



2; — 



27^ 2 



(D.24) 



In order to compute the asymptotic behavior of (D.24), we use the Laplace method. The 



maximum of the concave function (p(z) = 1 — z — (4/27) z 2 falls at c = 2/3, with 0(c) = 
and 4>"(c) = —9/2. Since, for any a < c < b, 



dzf(z)e 



T<t>(z) 



T<t>(c) 



(D.25) 
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the asymptotic behavior of / reads: 



2wy/3- 



(D.26) 



and 



(3AG{V) ~ -PgV + 6f3aV 2/3 + ln(6/3aF 2/3 ) - ln(27rV3) . 



(D.27) 



The last two terms in Eq. (D.27) give the subleading corrections to CNT as formulated for 



cubic clusters. The quality of the approximation (D.27) can be judged from the inset of 



Fig. 8, which shows a very good matching for all V^'s except for the smallest values, similarly 
to what occurs for rectangles (cf. Fig. 4). 

The calculation of AG can also be performed when a further energy cost, v per unit 



length, is assumed for the edges. Equation (D.16) is then modified to 



0AG{V) = -PgV + 6f3aV 2/3 + 120uV 



1/3 



In 



+oo 



dxdy exp {-^ (e x+y + e~ x + e~ y -3)-^ (e~ x - y + e x + e y - 3) ) 



(D.28) 



where n = 6/3aV 2/3 and r 2 = !2/3vV 1/3 . By the same line of reasoning as followed above 
we arrive at f3F s = — In I(V) with 



1 = 2 dy 



- x e -(n/3)v-75/(3v) / 2 Tl 2r 2 
K 



y 



(D.29) 



'o V V 3 v^ 3 

Laplace method can still be invoked to extract the asymptotic behavior of J, which turns 
out to be 

- In I ~ n + r 2 + ln(n + r 2 ) - ln(27rv / 3) • (D.30) 
From the above formula, we get 

/3AG(V) ~ + 6/3aV 2/3 + 12/3z/y 1/3 + \n(Qf3aV 2/3 + \2PvV l/3 ) - ln(27rV3) . (D.31) 



In Fig. 8, we compare the approximation (D.31) with the exact value. We see that the 
agreement is good for not too small V. 
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4. Ellipsoids 



Let us finally study the case of an ellipsoidal cluster. Volume and area of an ellipsoid 
with semiaxes a, b, and c are respectively given by 

4 /6c 2 \ 

V = -nabc and A = 2tt c 2 + FQ|m) + bVa 2 - c 2 E(<j>\m) , 

3 V v a, 2 — c 2 / 

(a>b>c and a > 6 > c ; A = 47rc 2 for a = 6 = c) (D.32) 

where 

m = T7w q ^ = i VT^ < 1 and = arcsin . (D.33) 

cr(er — cr) l — c z /a z a 

F and E are elliptic integrals of the first and second kind, respectively. For — 7r/2 < <fi < tt/2, 

they are defined as 

r<t> i r<t> 

F(<f)\m) = / dx — ^^^=^^^ and i£(0|m) = / dx y 1 — m sin 2 x . (D.34) 
</o vl-m sin 2 x ./ o 



Let now A(a, b, c) be the surface area of an ellipsoid of semiaxes a, b, and c (not necessarily 
in descending order). By the usual transformations, the surface free energy becomes 

" + °° -ex l-(3a(—X /3 

ab ^ ] V 47T / V ' ' ab J I 3 



./^(V) = -In // dad6^exp^-/3a( ^ A f a, 6 1 H + In 4 * 



o 



r 



In 



-co 



— oo 



dqdp exp J -— (A (e q , e p , e~ q - p ) - 4tt) \ + In — , (D.35) 

I 47T J 3 



where r = /3a(367r) 1/3 y 2/3 . To obtain pAG(V), it is sufficient to add -/3gV to ( |D.35fr . For 
(3g = (3<j = 1, the plot of this function is reported in Fig. 9. In the same figure, /3AG(V) is 
compared with the asymptotic estimate 

(3AG(V) ~ -pgV + /3(x(367r) 1/3 V 2/3 + ln(/3a(367r) 1/3 ^ 2/3 ) - 0.4849 , (D.36) 

where the last two terms give the correction to CNT as formulated now for spherical clusters. 
Judging from the inset of Fig. 9, which shows the difference between the approximate and 



exact values of /3AG(V), the estimate (D.36) is very good for all V's except for the very 
small ones. 



The strong similarity between (D.36) and (D.27), together with the high accuracy with 



which they reproduce the profile of (3AG(V) for ellipsoids and cuboids respectively, indicates 
that the difference between envisaging the nucleus as ellipsoidal rather than cuboidal entirely 
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FIG. 9: (Color online). Ellipsoids: /3AG(V) vs. V for pg = j3a = 1 (black solid line). The red 



dotted line is the approximant (D.36). Inset, the difference between the large-V estimate (D.36) 
and the exact @AG(V). 

lies in the value of er, which for an ellipsoid is 6/(367r) 1 / 3 ~ 1.241 times the cuboidal one. 
This occurs exactly as in CNT where the same relation holds between the values of a for 
spheres and cubes. 
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